import numpy as np
from collections import namedtuple
from scipy.optimize import fmin, brute, brent
from matplotlib import pyplot as plt
from pdb import set_trace
"""
    Module to calculate the drying dynamic in one dimension.
    The solution is based on the "phase change front" approach
    first proposed by Plank. 
    In this approach, the tempearture at the front is determined by
    equating the energy flux due to heat transfer to the energy flux
    due to solvent evaporation (plus the transient sensible heat
    increase of the bed). That is
    q_wall(T) = j(T)*H_vap + <cP>dT;
    where,
    q_wall(T) ~ Tf - T
    j(T) ~ (P_sat(T) - P) / (delta z)


                        |   ^ solvent flux  |
                        |   |               |
                        |-.-.-.-.-.-.-.-.-.-| <- bed height surface
       0-flux boundary-{|                   |
                        |                   |
                        |...................| <- solvent front (z meters)
                        |                   |
                        ---------------------
                            ^
                            |
                            heat flux
    References:
        Michaud et al Chemical engineering research and design 86 (2008) 606-611
        Kohout et al Mathematical modeling of solvent drying from a static particle
            bed Chemical engineering science (2006) 3674-3685

    Usage:
    simulate(parameters)
        parameters is a dictionary of physical parameters (see function simulate)
"""
            
def vapor_pressure(T, A, B, C):
    """
    Calculate the vapor pressure at temperature, T, using
    parameters for the Antoine eq. and return the vapor pressure in Pa.
    T - temperature to calculate the vaport pressure (K)
    parameters - Antoine parameters 
    return P_vap - vapor pressure at T (Pa) 
    """
    T -= 273
    P_vap = 10**(A-B/(C+T)) * 1.01325*10**5/760
    return P_vap

def mass_flux(z, P, z_c, P_v, K, h, r):
    """
    The vapor mass flux perpedicular to the front. The mass flux is equal to
    the vapor volumetric flux times the vapor density. The vapor volumetric
    flux is calculated from Darcy's Law.     
    volumetric flux = velocity = -(permeability/viscosity) * grad(pressure)
    m^3/(m^2*s) = m/s =  m^2 / Pa s *  Pa / m

    The pressure gradient is calculated as the difference between the vessel
    pressure and the pressure at the front divided by the distance from the
    front to the cake surface.
    grad(pressure) = (P_front - P_surface) / (z_front - z_surface)

    The mass flux is equal to the volumetric flux times the vapor density.
    mass_flux = density * velocity
    kg/(m^2*s) = kg/m^3 * m/s  

    P : pressure at z (Pa)
    z : axial distance from the cake bottome (m)
    P_v : vessel pressure (Pa)
    z_c : cake height (m)
    K : vapor phase permeability for Darcy's law (m^2)
    h : vapor phase viscosity (Pa s)
    r : vapor phase density (kg/m^3)
    return j : mass flux (kg / m^2 s)
    """
    j = - r * (K/h) * (P_v - P) / (z_c - z)
    return j

def Newton_cooling_parameter(X, h_w, h_f, X_cr1, X_cr2):
    """
    Calculate Newton's cooling parameter, h from the heat transfer
    fluid to the bed (i.e. jacket fluid -> wall -> bed).
    As the solvent concentration decreases the heat transfer
    coefficient from the wall to the bed decreases.
    Michaud Optimization of crystalline powders vacuum contact drying with
    intermittent stirring, Chemical engineeering research and design 86 (2008)
    606-611

    X : actual solvent concentration
    h_w : heat transfer coefficient from wall to cake (W / m^2 /K)
    h_f : heat transfer coefficient from heating fluid to wall (W / m^2 /K)
    X_cr1 : first critical solvent concentration
    X_cr2 : second critical solvent concentration
    return h : Newton cooling parameter (W / m^2 /K)
    """
    if X > X_cr1:
        h = h_f
    elif X < X_cr1 and X > X_cr2:
        h_g = 1 / (1/h_f + 1/h_w)
        h = (h_f - h_g)/(X_cr1 - X_cr2) * (X - X_cr2) + h_g
    else:
        h = h_w
    return h

def jacket_temperature(t, Tf_i, grad_Tf, Tf_s):
    """
    Calculate the jacket temperature for transient heat changer fluid
    temperature, to simulate the transient response of 
    starting the heat exchanger.
    t : time since the heat exchanger was started (typically start of experiment) (s)
    Tf_i : initial heat exchanger fluid temperture (K)
    grad_Tf : rate of change of heat fluid temperature (K / s)
    Tf_s : set point heat exchanger temperature (K)
    return : T_f the temperature of the heat fluid at time t (K)
    """
    T_f = np.min([Tf_i + grad_Tf*t,Tf_s])
    return T_f

def heat_flux(T, z, h, T_f, l_d):
    """
    Calculation of the energy flux through the wall into the filter cake. The
    flux is calculated as the effective heat transfer coefficient from the
    heating fluid through the wall through (any) dry bed to the front, divided
    by the difference between the temperature of the heat fluid and the
    temperature of the dry bed front.
    T : temperature of front (K)
    z : location of front (m) 
    T_f : temperature of fluid (K)
    h : Newton cooling parameter (W / m^2 /K)
    l_d : thermal conductivity dry layer (W / m /K)
    return q_w : heat flux (W / m^2)
    """
    q_w = (1 / (1/h + z/l_d)) * (T_f - T)
    return q_w

def mixture_heat_capacity(X, cp_s, r_s, cp_l, r_l, Y):
    """
    Calculate the volume average heat capacity of the solid/solvent mixture per
    unit area of vessel cross section. This calculations assumes that the heat
    capacity of the gas in the void volume of the solids is zero.
    X : mass fraction of solvent
    cp_s : heat capacity solid (J/kg)
    r_s : density solid (kg/m^3)
    cp_l : heat capacity liquid (J/kg)
    r_l : density liquid (kg/m^3)
    return heat_capacity : heat capacity of solid/solvent mixture (J/kg)
    """
    heat_capacity = cp_l * r_l * (X / (X+1)) + cp_s * r_s *(1 / (X+1))
    return heat_capacity

def vapor_density(T, P_v, mw):
    """
    Calculate the vapor density from the ideal gas law.
    T : Temperature (K)
    P_v : local pressure (vessel pressure) (Pa)
    mw : molecular weight (g/mol)
    return vapor density (kg/m^3)
    """
    r = (P_v * mw ) / (8.314 * T)  
    return r

def local_temp_objective_function(T, p, t, T_i, z, X, h, T_f): 
    """
    Objective funtion to be minimized to solve for local temperature.
    The conductive heat flux 
    heat_flux(T) = mass_flux(P_sat(T))*H_vap + rate change_sensible_heat(T_i-T)
    [(mass_flux*H_vap + change_sensible_heat(T_i-T) - heat_flux)(T)]^2 = 0
    T : temperature to solve for steady state (K)
    p : namedtuple of parameters used in the simulation
    z : front position (m)
    r : vapor density (kg/m^3)
    """
    q = heat_flux(T, z, h, T_f, p.l_d)
    P_sat = vapor_pressure(T, p.A1, p.A2, p.A3)
    r = vapor_density(T, p.P_v, p.mw)
    j = mass_flux(z, P_sat, p.z_c, p.P_v, p.K, p.v, r)
    """ Heat capacity of the wet portion of the bed """
    #S_w = mixture_heat_capacity(X, p.cp_s, p.r_s, p.cp_l, p.r_l, p.Y) * (p.z_c - z) 
    """ Heat capacity of the dry portion of the bed """
    #S_d = mixture_heat_capacity(0, p.cp_s, p.r_s, p.cp_l, p.r_l, p.Y) * z
    """ The total power (energy per unit time). Note the 0.5 for the
    dry is due to the gradient: <T> ~ 1/2(T_wall - T_front) """
    #S_t = (S_w * (T - T_i) + S_d * 0.5 * (T - T_i) ) / p.dt
    """ Energy balance """
    #objective_function = (j * p.H_vap + S_t - q)**2
    objective_function = (j * p.H_vap - q)**2
    return objective_function

def simulate(parameters):
    """
    Psuedo steady state Agitated Filter Drying model in one dimension
    See Kohout et al Chem. Eng. Sci. 61 (2006)
    Developed by: Brent Maranzano
    Created date: 18, May 2011
    Last modified date: 27, June 2011
    params : dictionary with following key
        r_s : solid density (kg/m^3)
        Y : dry solid porosity (V_void)/V ; where V_void=V-V_solid)
        l_d : thermal conductivity dry layer (W/m/K)
        cp_s : solid heat capacity (J/kg K)
        r_l : liquid density (kg / m^3)
        cp_l : liquid heat capacity (J / kg K)
        v : vapor viscosity (Pa s)
        mw : vapor molecular weight (g/mol) 
        H_vap : enthalpy of vaprization (J/mol) #TODO check units
        A1, A2, A3 : Antione equation parameters to give vapor pressure in mmHg
        K : vapor permeability (m^2)
        h_w : outer wall newton cooling parameter (W/m^2/K)
        h_f : inner wall newton cooling parameter (W/m^2/K)
        X_cr1 : first critical solvent concentration 
        X_cr2 : second critical solvent concentration 
        z_c : cake height (m)
        D : vessel diameter (m)
        P_v : vessel pressure (Pa)
        Tf_i : initial temperature heat exchange fluid (K) 
        grad_Tf : slope change heat exchanger fluid (K/s) 
        Tf_s : set point of the heat exchanger fluid (K) 
        To : initial temperature (CK 
        Xo : initial mass fraction
        dt : integration time increment (s)
    """
    for k, v in parameters.iteritems():
        parameters[k] = float(v)

    """ create a named tuple of the parameters, so can use dot notation """
    ParamStruct = namedtuple('ParamStruc', parameters.keys())
    p = ParamStruct(*tuple(parameters.values()))

    """ Define the cross section area of vessel (A[=]m^2) and the total mass of
    solid (m_s[=]kg). """
    A = np.pi * p.D**2 / 4.
    m_s = (1-p.Y) * A * p.z_c * p.r_s

    """ set intitial time (t=0 s), evaporation rate (m=0 kg/s), front position
    (z=0 m), and temperature (T=parameters['To']) """
    t = 0                                 
    m = 0                                 
    z = 0                                 
    X = p.Xo #p.Y / (1-p.Y) * p.r_l / p.r_s     
    T = p.To

    results = [dict(t=t, m=m, z=z, X=X, T=T)]
    """ begin calculation """
    while X > 0:
        """ Calculate some params outside objective function for speed improvement. """
        h = Newton_cooling_parameter(X, p.h_w, p.h_f, p.X_cr1, p.X_cr2)
        T_f = jacket_temperature(t, p.Tf_i, p.grad_Tf, p.Tf_s)
        try:
            """
            T, fopt, iter, funcalls, warnflag = fmin(local_temp_objective_function, 
                   T, (p, t, T, z, X),
                   xtol=0.0001, ftol=0.0001, maxiter=10000, maxfun=None, 
                   full_output=1, disp=0, retall=0)
            opt = brute(local_temp_objective_function, (slice(280,p.Tf_s),),
               args=(p, t, T, z, X), full_output=True, finish=fmin)
            """
            opt = brent(local_temp_objective_function, 
                    args=(p, t, T, z, X, h, T_f), 
                    brack=(280,p.Tf_s), full_output=True)
        except KeyboardInterrupt:
            return results

        T = opt[0]
        P_sat = vapor_pressure(T, p.A1, p.A2, p.A3)
        r = vapor_density(T, p.P_v, p.mw)
        j = mass_flux(z, P_sat, p.z_c, p.P_v, p.K, p.v, r)
        m = -1 * j * A
        X += (m/m_s) * p.dt 

        if X < p.X_cr2:
            """ the front begins to move """
            z +=  j / (p.X_cr2 * (1 - p.Y) * p.r_s) * p.dt

        results.append({'t':t, 'm':m, 'z':z, 'X':X, 'T':T})
        t += p.dt

    return results

def graphical_check(parameters, t, z, X):
    """ 
    For the given physical properties (parameters), current front position,
    and solvent concentration, make a plot of the thermal heat flux (energy in)
    as a function of temperature and the mass flux * enthalpy of vaporizaion
    (energy out) as a function of temperature on the left axis.  The two curves
    should intersect near the steady state value (differenced by the sensible
    heat increase).  Also plot the vapor pressure as a function of temperature
    on the right axis. The plots aid in trouble shooting poor convergence results.

    parameters : dictionary of physical properties 
    t : simulation elapsed time (s) for calculating heat fluid temperature
    z : current front position (m)
    X : current solvent concentration
    """
    ParamStruct = namedtuple('ParamStruc', parameters.keys())
    p = ParamStruct(*tuple(parameters.values()))

    results = np.zeros((200), dtype={
        'names':['T', 'q', 'j', 'Psat'],
        'formats':['float32'] * 4})

    T_f = jacket_temperature(t, p.Tf_i, p.grad_Tf, p.Tf_s)
    h = Newton_cooling_parameter(X, p.h_w, p.h_f, p.X_cr1, p.X_cr2)

    for i, T in enumerate(np.linspace(20, 350, 200)):
        q = heat_flux(T, z, h, T_f, p.l_d)
        P_sat =vapor_pressure(T, p.A1, p.A2, p.A3) 
        r = vapor_density(T, p.P_v, p.mw)
        j = mass_flux(z, P_sat, p.z_c, p.P_v, p.K, p.v, r) * p.H_vap
        results[i] = T, q, j, P_sat

    fig = plt.figure()
    axl = fig.add_subplot(111)
    lines = axl.plot(results['T'], results['q'], color='red', label='q')
    lines = axl.plot(results['T'], results['j'], color='blue', label='j*H_vap')
    axr = axl.twinx()
    lines = axr.plot(results['T'], results['Psat']*1.012e-5, color='green', label='Psat')
    axl.set_xlabel('T (K)')
    axl.set_ylabel('q (W), j*Hvap (W)')
    axr.set_ylabel('Psat (bar)')
    axl.legend()
    plt.draw(); plt.show()
    return results

if __name__=='__main__':
    parameters = dict(r_s=2100., Y=0.45, l_d=0.169, cp_s=837.,
        r_l=985., cp_l=4186.,
        v=8.6*10**-6, mw=30., H_vap=9.186*10**5,
        A1=8.20417, A2=1642.89, A3=230.3,
        K=0.8*10**-12, h_w=10., h_f=10.,
        X_cr1=0.1, X_cr2=0.08,
        z_c=0.0065, D=0.003, P_v=7000., Tf_i=298., grad_Tf=0.0002, Tf_s=318.,
        To=293., Xo=0.2, dt=1.)
    simulate(parameters)
